#!/usr/bin/perl 
#===============================================================================
#
#         FILE: MatchPeak.pl
#
#        USAGE: ./MatchPeak.pl  
#
#  DESCRIPTION: The script takes output from tandem_ms_ 
#
#      OPTIONS: ---
# REQUIREMENTS: ---
#         BUGS: ---
#        NOTES: ---
#       AUTHOR: Han Hu (hh), hh1985@bu.edu
#      COMPANY: Bioinformatics Program & Center for Biomedical Mass Spectrometry, Boston University
#      VERSION: 1.0
#      CREATED: 8/18/2011 9:49:43 AM
#     REVISION: ---
#==============================================================================

use strict;
use warnings;
use Getopt::Std;

my %Options;
getopt('pxis',\%Options);
my ($pre, $ppm);
my $dirname = "./data";
my $libdir = "./output";
if(!defined($Options{x})) {
	die "$pre is not provided!\n";
} else {
	$pre = $Options{x};
}
if(!defined($Options{p})) {
	$ppm = 5e-6;
} else {
	$ppm = $Options{p};
}
my $outfile = "matched_" . $pre . ".txt";
my $infile = "data_" . $pre . ".dta";
my $libfile = "gag_library_$pre.csv";

open(LIB, "$libdir/$libfile") or die "Cannot read the file:$!\n";
my (@library, @reconstruct);
my $line_id = 0;
while(my $line = <LIB>) {
	
	chomp $line;
	my @group = split /\t/, $line;
	push @library, \@group;
	my @all_loss;
	my @neutral_loss = split "-", $group[7];

	if(defined($group[8])){
		my @internal_loss = split "-", $group[8];
		## Also consider about the situation with no loss.
		@all_loss = (0, @neutral_loss, @internal_loss);
	} else {
		@all_loss = (0, @neutral_loss);
	}

	foreach my $loss (@all_loss){
		## Create the track link from mass to line number.

		#push @reconstruct, [($group[2]-$loss), $line_id];
		if($group[2] < $loss){
			next;
		}
		push @reconstruct, [$line_id, $group[2], $group[2]-$loss, $loss];

	}	
	
	$line_id++;
}
close(LIB);

open(PK, "$dirname/$infile") or die "Cannot read the file:$!\n";
chomp(my $top = <PK>);
my @pre = split /\t/, $top;
my @pro;
while(my $line = <PK>) {
	chomp $line;
	my @group = split /\t/, $line;
	next if $group[2]==0; ## If the charge state is 0, it cannot be identified as true.
	push @pro, \@group;
}
close(PK);

## Match mass against the library.
my @mass = map {$_->[2]} @reconstruct;
## Sort numerically ascending 
my @sorted_lib = sort {$a <=> $b} @mass;
## Get the corresponding order 
my @order_lib = sort {$mass[$a] <=> $mass[$b]} 0..$#mass;

my @peak_mass = map {$_->[2]} @pro;
my @sorted_peak = sort {$a <=> $b} @peak_mass;

use SequenceToFragments;

my ($pt_lib,$pt_data) = (0,0);
my @hit_mass;

open(FL, ">$dirname/$outfile") or die "Cannot write to the file $outfile:$!\n";
foreach my $pms (@sorted_peak) {
	## Quickly locate the position.
	for(my $j=$pt_lib; $j<@sorted_lib; $j++){
		my $error = abs($pms-$sorted_lib[$j])/$sorted_lib[$j];
		## Consider about mass shift.
		## Hit one!
		if($error <= 2 * $ppm){
			my $LineNum = $reconstruct[$order_lib[$j]]->[0];
			my $string = join "\t", ($pms, @{$library[$LineNum]}[0..6], @{$reconstruct[$order_lib[$j]]}[2,3]);
			print FL "$string\n";	
		} elsif($sorted_lib[$j] > $pms) {
			## Imporsible to get the hit.
			$pt_lib = $j;
			last;
		}
	}
}

close(FL);

